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Abstract One of the main purposes for conducting well tests is to obtain information abont reser¬ 
voir parameters, such as its permeability and the existence of skin effects. Determining individual 
layer properties from pressure transient data in multilayer reservoirs is challenging, since pressure 
behavior is influenced by the properties of all layers. This work extends two existing methods for 
estimating individual layer properties in multilayer systems with vertical wells to stratified reser¬ 
voirs with multilateral horizontal wells. These methods are the rate-normalized pressnre analysis 
and the delta transient method. Both methods reqnire a clear identification of radial flow regimes, 
which may not occur in a practical case. We show that the Nelder-Mead global method algorithm 
can also be used to evaluate layer permeabilities and skin. This method represents an alternative to 
estimate layer parameters in cases where radial flow regimes are difficulty to identify or nonexistent, 
since it does not depend on the occurrence of any specific flow regime. The proposed techniques 
are applied on a pair of synthetic cases, where pressure and layer flow-rate profiles are compnted 
from a previously existing analytical model. Resnlts show that all three methods are able to yield 
good estimates for layer horizontal and vertical permeabilities, but only the Nelder-Mead method 
is able to accurately estimate the individual skin factors. Moreover, the Nelder-Mead method can 
provide good estimates for layer properties even if late-time radial flow is not observed. Nonethe¬ 
less, our results suggest that early-time data is critical to estimate layer vertical permeabilities and 
mechanical skin. 

Keywords Well Testing • Multilateral Wells • Stratified Reservoirs • Rate-normalized pressure 
analysis • Delta transient method • Nelder-Mead algorithm 


1 Introduction 

Reservoir characterization is a major goal of well testing. Based on the pressnre response measured 
during the test, it is possible to estimate reservoir features such as permeability, average pressure, 
skin factor and outer boundary condition. In multilayer systems, it is desirable to obtain information 
regarding each individual layer, so that the reservoir productivity may be more accurately estimated. 

In multilayer stratified systems, conventional interpretation techniques allow the determination 
of the reservoir eqnivalent permeability and skin (Cobb et ah, 1972; Raghavan et ah, 1974; Larsen, 
1982; Kuchuk et ah, 1986a; Ehlig-Economides and Joseph, 1987; Raghavan, 1989; Kuchnk and 
Wilkinson, 1991). Computing individual layer properties, however, is not so straightforward, since 
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pressure response at the wellbore is simultaneously influenced by all layers (Lefkovits et ah, 1961; 
Cobb et al., 1972; Kuchuk et ah, 1986b; Guimaraes and Galvao, 2017). 

Larsen (1982) applied the analytical approximation proposed by Larsen (1981) for pressure 
behavior in two-layered commingled reservoirs with vertical wells to obtain estimates for layer flow 
capacities and skin. Layer properties are determined through a trial-and-error procedure that aims to 
match the observed pressure data and the theoretical pressure response obtained using the analytical 
expression developed by Larsen (1981). If data are sufficiently smooth, then layer properties may 
be evaluated using analytical expressions, instead of a trial-and-error procedure. Larsen (1982) also 
stated that this procedure could be applied in reservoirs with crossflow, provided that data after 
crossflow effects become relevant are avoided. 

A method to evaluate layer properties in reservoirs with unequal initial pressure was derived 
by Aly et al. (1994). Their method, which is called Derivative Extreme Method, is based on the 
pre-production pressure data, rather than conventional well testing. It consists of identifying the 
point where pressure derivative attains its extreme (maximum or minimum) value. Then, a system 
of equations is built so that layer properties such as permeabilities, skin factors and porosities may 
be computed. 

Gao (1987) proposed a means to determine layer permeabilities in a two-layer reservoir with 
vertical wellbore from well testing data. The method consists of isolating the layers with packers so 
that individual well test analysis may be performed for each layer. Thereby, conventional single-layer 
interpretation methods may be applied to obtain layer permeabilities and skin factors. 

As pointed out by Raghavan (1989), isolating the layers and running individual pressure transient 
tests for each one of them is often unfeasible in practice. Therefore, Raghavan (1989) claimed that 
the rate-normalized pressure analysis (or RNPA) represents a more viable alternative to determine 
individual layer properties. This technique was first introduced by Ehlig-Economides and Joseph 
(1987), in a work that also provided a detailed description of the analytical model for single-phase 
flow in multilayer reservoirs with vertical wells. RNPA, however, presents an intrinsic operational 
setback, as it requires several measurements of layer flow-rates. 

An alternative solution was developed by Guimaraes and Galvao (2017). In their work, a new 
layer parameter denoted as ’’delta transient” is defined. Based on production logging data from 
merely two distinct time steps, the delta transient parameter is determined and, hence, so are layer 
permeabilities and skin factors. Their work presents a strong operational drive by requiring only 
a few hours of rig time to acquire the test data. In a further work, Galvao and Guimaraes (2017) 
extended the delta transient method (or DTM) to tests with a multiple flow-rates scheme. 

All the methods mentioned above present one common feature; they were designed for multilayer 
reservoirs with vertical wells. Thereby, current approaches to estimate individual layer parameters 
are not suitable for reservoirs with horizontal wells. Interest in drilling horizontal wells has been 
increasing, as they present some benefits compared to vertical wells. They mitigate fluid coning 
effects and their performance is superior in reservoirs with high vertical permeability (Clouts and 
Ramey Jr., 1986; Ozkan et ah, 1989; Algharaib et ah, 2006). Moreover, one horizontal well may pro¬ 
duce the same amount of oil as several vertical wells (Joshi, 2003). Horizontal wells also present high 
productivity, due to their larger contact area with the reservoir (Algharaib et ah, 2006; Biryukov and 
Kuchuk, 2015). However, to the best of our knowledge, there is currently no reliable interpretation 
technique to obtain estimates for individual layer properties from pressure data collected during 
well testing in stratified reservoirs with multilateral horizontal wells. 

Thus, the main purpose of this work is to achieve a means of determining individual layer 
permeabilities and skin factors from well testing data in multilayer reservoirs with horizontal well¬ 
bores. In this paper, the RNPA (Ehlig-Economides and Joseph, 1987; Raghavan, 1989) and the 
DTM (Galvao and Guimaraes, 2017; Guimaraes and Galvao, 2017) are extended so that they can 
be applied in horizontal wells also. Both RNPA and DTM demand nothing more than a general 
worksheet software to be performed. 

The proposed extensions of RNPA and DTM rely on the identification of the two distinct radial 
flow regimes that might occur during single-phase flow in horizontal wells: one at early times, while 
the vertical boundaries are not felt by the wellbore and another at late times, when pressure diffusion 
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propagates radially in the horizontal plane (Rosa and Carvalho, 1989; Jelmert and Thompson, 1991; 
Knchuk, 1995; Nie et ah, 2011). 

Yet, depending on the reservoir properties and test duration, sometimes no radial flow regime 
may be observed. In those cases, optimization methods consist of an alternative to pressure transient 
interpretation techniques. Reservoir parameters are estimated by minimizing an objective function, 
so that the computed parameters yield pressure and flow-rate responses that match the observed 
data. In this context, gradient-free methods are usually desired, since they are more computationally 
efficient. The Nelder-Mead method (also referred to as simplex or polytope method) is a global 
optimization algorithm that aims to find a set of parameters that minimize the objective function 
by sequentially applying four simple operations: reflection, expansion, contraction and shrinkage 
(Nelder and Mead, 1965; Barati, 2011; Feng et ah, 2020). It may be used to estimate reservoir 
parameters by matching pressure data (Rahmati et ah, 2013) when coupled with analytical models. 

There exist analytical solutions for pressure behavior in single-layer reservoirs with horizontal 
wells both in Laplace domain (Goode and Thambynayagam, 1987; Kuchuk et ah, 1991) and in real 
space (Daviau et ah, 1988; Ozkan et ah, 1989). Analytical models for multilateral wells in multilayer 
reservoirs with multilateral horizontal wells are also available (Vo and Madden, 1995; Larsen, 2000; 
Yildiz, 2003; Pan et ah, 2010). Layer flow-rate profiles may also be determined from those models. 
Thus, a third approach was conducted by using the analytical solution from Pan et al. (2010) as 
forward model for the Nelder-Mead global optimization algorithm (Nelder and Mead, 1965). 

The paper is organized as follows: Section 2 presents the early- and late-time logarithmic ap¬ 
proximations for pressure change, whereas Section 3 depicts the proposed techniques to determine 
layer permeabilities. The obtained results are displayed in Section 4. Finally, Section 5 presents the 
main conclusions of this work. 


2 Logarithmic approximations for pressure change 

The real space solution for single-phase flow of a slightly compressible fluid with constant com¬ 
pressibility and viscosity in reservoirs with multilateral horizontal wells was developed by Pan et al. 
(2010). Their model is outlined in Appendix A. Figure 1 displays a schematic of the reservoir model 
considered in this study. Eq. (A-1) and the linear system of equations given by Eq. (A-4) (given 
in Appendix A) apply at all times. Nonetheless, early- and late-time approximations can be de¬ 
rived, considering the distinct flow regimes observed in horizontal wells. Here and throughout, all 
equations are expressed in terms of consistent unit systems. 



Fig. 1: Schematic of the Reservoir Model Considered in This Study 
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As evidenced in Fig. 1, layer properties such as permeability and thickness may be different in 
each layer. The considered model also accounts for distinct wellbore lengths and well off-centering. 
Therefore, the smallest distance between the wellbore and a vertical boundary (denoted as dzj) 
may be measured from the wellbore to the layer top boundary or to the layer bottom boundary, 
depending on the wellbore position. 


2.1 Early-time approximation 


Immediately after the wellbore is open, equipotential curves are concentric ellipses in the vertical 
plane perpendicular to the wellbore. A radial flow is observed, and pressure response is similar to 
a vertical well in an infinite acting reservoir whose thickness equals the horizontal wellbore length 
(Daviau et ah, 1988; Rosa and Carvalho, 1989; Nie et ah, 2011). While early-time radial flow occurs, 
an approximation for pressure change at the wellbore in a multilayer reservoir with a multilateral 
horizontal well is given by: 


^Pwf (1) 


qBjj, 

Lf {^zx')eq 


'1 / 4(fc^a;)egt \ 

2 Vexp(7)li(<(>Ct) 


+ 5e 


( 1 ) 


where q stands for the production flow-rate, B stands for formation volume factor, /i stands for 
oil viscosity, Lt stands for the wellbore total length and r-w stands for the wellbore radius. The 
equivalent permeability in the vertical plane, denoted as (kzx)eq, is defined as: 


X) Ljkzxj 

{kzx)eq = - z , ( 2 ) 

Lit 

where kzxj = \/kxyjkzj is the geometric mean permeability in the 0 a;-plane for layer j. Besides, in 
Eq. (1), {<j)Ct)eq represents a thickness weighted average of the porosity-compressibility product: 


n 

X) 4>jCtjhj 

{4>Ct)eq = - • (3) 

In Eq. (1), Seq denotes the equivalent mechanical skin factor. Analogously to the definition 
proposed by Ehlig-Economides and Joseph (1987), it may be computed as a weighted average of 
the skin factors in each layer: 


n 

X) kzxjLjSj 

“ {kzx)eqLt 

Eq. (1) is valid as long as the early-time radial flow is simultaneously occurring in all layers. 
Pressure change during early-time radial flow increases linearly with the natural logarithm of time. 
Thus, the equivalent permeability in the vertical plane {kzx)eq can be estimated by taking the 
pressure derivative with respect to the natural logarithm of time (hereafter denoted simply as 
pressure derivative) (Yildiz, 2003): 
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( 5 ) 


Thereby, after {kzx)eq has been determined from Eq. (5), the equivalent mechanical skin may 
be computed as: 


Sea = X 


Apixf (terf ) 


^erf 


- \n{terf) - In 


4i(kzx}e 
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( 6 ) 


where terf is any chosen time during the early-time radial flow and Ap^ filer f) is its corresponding 
pressure change. 
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2.2 Late-time approximation 


At late times, a pseudo-radial (or late-time radial) flow develops, as flow in the vertical direction 
becomes negligible and pressure diffusion propagates predominantly in the horizontal plane (Rosa 
and Carvalho, 1989; Odeh and Babu, 1990; Kuchuk et ah, 1991; Nie et ah, 2011). Once the late-time 
radial flow prevails in all layers, pressure response at the wellbore may be approximated by: 


f (^) 


( ^{kxy)eqt \ „ 

ht{kxy)eq [2 \exp{-/) eqrl J ’ 


( 7 ) 


where ht stands for the reservoir total thickness. The equivalent permeability in the horizontal plane 
{kxy)eq is defined as: 


X) kjkxyj 

{kxy)eq = “ . (8) 

At late times, besides the mechanical damage, there is also a pseudo-skin effect due to the 
convergence of stream lines beyond the well tips (Goode and Thambynayagam, 1987; Kuchuk et ah, 
1991). Therefore, in Eq. (7), Steq stands for the equivalent total skin factor. Analogously to Eq. (4), 
it consists of a weighted average of total skin factors in all layers: 


n 

X] kxyjhjStj 

Q _ _ 

(kxy)eqht • 

where Stj stands for the total skin factor in layer j and is defined as (Daviau et ah, 1988): 


( 9 ) 


c — c I c 
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( 10 ) 


where Sgj stands for the pseudo-skin in layer j. For a laterally infinite reservoir, it is computed as 
(Daviau et ah, 1988): 


^gj — In 



kxyj kj 
kzxj Lj 


( 11 ) 


Eq. (7) shows that pressure change is a linear function of the natural logarithm of time. Therefore, 
pressure derivative during late-time radial flow may be used to evaluate the equivalent permeability 
in the cry-plane (Yildiz, 2003): 


dApxuf _ _ qBp 

d\n{t) 2ht{kxy)eq 
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qBp. 
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( 12 ) 


After the horizontal permeability has been computed, the equivalent total skin factor may be 
determined as: 


St 


eq 


1 r ^Pu,f{tlrf) _ _ /- A{kxy)eq - \ 

2 L mirf \exp(-y)fi((j)Ct)eqriJ 


(13) 


where tirf is an arbitrary time during the late-time radial flow and Ap-u,f{tirf) is its corresponding 
pressure change. 

The equivalent total skin factor defined in Eq. (9) represents a weighted average of the total 
skin factors in all layers, encompassing both stream lines geometry and mechanical damage around 
the well. However, since wellbore lengths and layer thicknesses may be different in each layer, the 
geometrical pseudo-skin is also different in each layer. Therefore, the use of late-time data to 
determine the mechanical damage from the equivalent total skin factor is quite challenging, and out 
of the scope of this work. 
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3 Proposed Methods for Estimating Layer Properties 

Eqs. (1) and (7) may also be applied for a given layer j. Thus, at early-times, pressure change in 
layer j may be approximated by: 


^Pwf (t) 


^zxj Lj 


-Inf-^ 


(14) 


The late-time logarithm approximation for pressure change in layer in its turn, is given by: 


^Pwf (^) 


^xyj hj 



4A)x'U7 i 
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To obtain Eqs. (14) and (15), it was assumed that wellbore pressure is the same in all layers, 
except for the hydrostatic column. Furthermore, it was considered that the durations of early- 
and late-time radial flow have been previously identihed from pressure derivative data. It is also 
important to highlight here that, in Eqs. (14) and (15), layer flow-rate qj is a function of time. As 
explained in Appendix A, flow-rates are updated at each time step and, therefore, may change in 
time (Pan et ah, 2010). 

Based on Eqs. (14) and (15), the RNPA (Ehlig-Economides and Joseph, 1987; Raghavan, 1989) 
and DTM (Galvao and Guimaraes, 2017; Guimaraes and Galvao, 2017) were extended to multilayer 
systems with multilateral horizontal wells, as the schematic portrayed in Fig. 1. 


3.1 Rate-Normalized Pressure Analysis 

Rate-normalized pressure change is dehned as the ratio between wellbore pressure and the instant 
flow-rate at a given time (Ehlig-Economides and Joseph, 1987; Raghavan, 1989). Starting from Eq. 
(14), the rate-normalized pressure change in layer j during early-time radial flow may be computed 
as: 


PNj {t) 


^Pwf (1) 


Bp-o 

^zxj Bj 


“^kzxjt 

2 



(16) 


The rate-normalized pressure analysis (RNPA) encompasses the changes in pressure and layer 
flow-rate into a single function of time. Therefore, it is more suitable for estimating individual layer 
permeabilities than merely the pressure data. Eq. (16) shows that the rate-normalized pressure 
change behaves as a linear function of the natural logarithm of time. Taking the derivative of Eq. 
(16) with respect to the natural logarithm of time: 


dpNjjt) _ _ Bfij _^ . _ Buj 

^ln(t) 2kzxjLj 2merfjLj 


(17) 


Thus, layer permeabilities in the 2 a;-plane may be obtained from Eq. (17), using the rate- 
normalized pressure data corresponding to the early-time radial flow. The mechanical skin factor in 
layer j, then, may be estimated as: 


S, = - 


PNj jterf) 
TYX^yf 


- ln(ter/) - In 


4fc, 


exp(7)<)>jMlCtjr2 


(18) 


where ter/ and pj^jiterf) represent, respectively, an arbitrary time point during which early-time 
radial flow is observed and the corresponding rate normalized pressure change in layer j. 

At late-times, an expression for the rate-normalized pressure change may be obtained from Eq. 
(15): 


PNj (t) 
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(19) 



Estimating Layer Properties in Reservoirs With Multilateral Wells 


7 


The horizontal permeability in layer j may be estimated from the derivative of Eq. (19) with 
respect to the natural logarithm of time: 


( 20 ) 


dpNjjt) _ _ Bpj _^ . _ Bpj 

c)ln(t) 2kxyjhj 2mirfjhj' 

Thus, for a given time point tirf such that late-time radial flow occurs, and its corresponding 
rate normalized pressure change PNj{tirf), total skin factor in layer j may also be computed by: 


Stj = I 


PNjjtlrf) 

TTlirf 


- - In 


Akx 


exp(7)0j^ctjr2 


( 21 ) 


3.2 Delta Transient Method 

The application of the RNPA presents an intrinsic operational drawback, which is the need to 
continuously measure layer flow-rates. To overcome this issue, Galvao and Guimaraes (2017) and 
Guimaraes and Galvao (2017) derived an alternative proposal to estimate individual layer perme¬ 
abilities and skin factors. They referred to their method as the Delta Transient Method (DTM). 
Their motivation to develop this procedure was to avoid the need of registering layer flow-rates sev¬ 
eral times. This section presents an extension of the DTM to multilayer reservoirs with multilateral 
horizontal wells. 

Based on the early-time logarithmic approximation for pressure change in layer j, given by Eq. 
(14), two variables 5erfj and /3j{t) shall be defined as follows: 

6,rfj = Hkzxj) + 25,; /?,(t) = In . (22) 

Thus, replacing Eq. (22) in Eq. (14): 


it) = [/?. (t) + ■ (23) 

ZiKzxj J-'j 

Rearranging Eq. (23) in a more convenient way: 

2kzxjLj _ Qjjt) if^jjt) + ^erfj] / 2 ^\ 

Bpj ^Pwfit) 

One should note that no further assumption was made to obtain Eq. (24), expect for those 
already mentioned in Section 2 and Appendix A. It requires only the logarithmic approximation 
given by Eq. (14) to be valid. Thus, for two distinct time points and t 2 during which early-time 
radial flow regime is occurring: 


er fj] 

Apxufitl) 


which implies that: 
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(26) 


(27) 


Thus, Serfj may be determined from Eq. (27). Then, the vertical permeability in layer j is 
estimated from Eq. (23): 
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2LjAp^f{t) 

Moreover, the mechanical skin Sj is evaluated from the definition of Serfj- 

Q _ ^erfj lll(^zxj) 


(28) 


(29) 


An analogous procednre may be applied nsing late radial flow data. The definition of 6, now, 
shonld be given by: 


^Irfj — ^^(^icyj ) T ‘2‘Stj~ 


(30) 


Then, similar computations show that the total skin factors and horizontal permeabilities in 
each layer may be obtained from late-time data: 


( <lipj ) 


( <ljPi ) 



t = t4^ 
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t = t3 

\ ^P^S ) 

t = t4 


(31) 


In Eq. (31), ts and t 4 denote two distinct time points after late-time radial flow has already 
started. Thus, horizontal permeability in layer j is determined as follows: 


2hjAp^f{ty 




Total skin factor is evaluated from the dehnition of Sirfj'- 


Su = 


_ ^Irfj lu(A:xyj) 


(32) 


(33) 


3.3 Reasonableness Check 

Galvao and Guimaraes (2017) claimed that a reasonableness check must be performed to assess if 
the computed layer properties are consistent with the estimated reservoir eqnivalent parameters. 
The reasonableness check consists of comparing the reservoir equivalent properties obtained from 
conventional drawdown/bnildup analysis with the compnted properties for each individual layer. If 
the difference between these results is relevant, then the assnmptions made in the development of 
RNPA and DTM (snch as the logarithmic approximations for pressnre change) are not valid (Galvao 
and Gnimaraes, 2017; Guimaraes and Galvao, 2017). In that case, the estimates for individual layer 
parameters should be discarded. 

In this work, we propose that the valnes for layer permeabilities evaluated from both the RNPA 
and the DTM shonld be applied in Eqs. (2) and (8) to compnte the reservoir equivalent permeabil¬ 
ities. These equivalent permeabilities shonld be compared to the reservoir eqnivalent permeabilities 
determined from Eqs. (5) and (12). 

Besides, the mechanical skin factors determined from early-time data should be applied in Eq. 
(4) and compared to the eqnivalent mechanical skin obtained from Eq. (6). Conventional pressure 
interpretation methods do not allow the determination of an equivalent mechanical skin based on 
late-time radial flow data, due to the pseudo-skin factor in each layer. Thereby, the reasonableness 
check for late-time skin factors should be performed using total skins, instead of the mechanical 
skin. Thns, the obtained estimates for Stj should be applied in Eq. (9) and compared to the value 
of Steq computed from Eq. (13). 

As stated by Guimaraes and Galvao (2017), if the equivalent properties compnted from the 
estimated individual layer properties are signihcantly different from those obtained from the con¬ 
ventional analysis, then the results should be discarded. In their work, Guimaraes and Galvao (2017) 
did not provide a more objective criteria to assess the resnlt of the reasonableness check. In this 
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work, the maximum acceptable difference between the equivalent parameters was set as 10%. This 
means that, if the estimated layer properties yield equivalent permeabilities and skin factors that 
are 10% higher (or lower) than the values obtained from conventional analysis, then the results 
should not be considered. 


3.4 Nelder-Mead Simplex Algorithm 


Application of the RNPA and DTM requires the constant derivative levels related to early- and late¬ 
time radial flow to be clearly identified. Otherwise, the logarithmic approximations given by Eqs. 
(14) and (15) are not valid. However, depending on layer thicknesses and wellbore laterals lengths, 
early-time radial flow may end so quickly that no constant derivative level is identified (Odeh and 
Babu, 1990; Kuchuk et ah, 1991). Identification of early-time radial flow is also compromised if the 
wellbore storage effect is relevant. Besides, if the wellbore is significantly long, late-time radial flow 
may not be observed during the test (Ozkan et ah, 1989; Rosa and Carvalho, 1989). 

Therefore, optimization methods represent an interesting alternative to estimate layer param¬ 
eters, since they do not rely on the identification of any flow regime. The Nelder-Mead simplex 
algorithm (NM) is a gradient-free global optimization procedure that is based on simply evaluating 
the objective function at the vertices of a simplex (Nelder and Mead, 1965; Feng et ah, 2020). The 
simplex is defined as a set of A" -|- 1 vertices, where N is the number of independent variables of the 
objective function. Then, the vertex that yields the highest objective function value is successively 
replaced until a convergence criteria is satisfied (Luersen and Le Riche, 2004; Barati, 2011). 

At each iteration, there are four possible operations that determine how the simplex should 
be updated: reflection, expansion, contraction and shrinkage. Then, after the proper operation is 
performed, a new simplex is generated and convergence is checked. If the convergence criteria is 
satisfied, the algorithm stops. Otherwise, a new iteration is performed (Ghiasi et ah, 2007; Rahmati 
et ah, 2013). A more detailed flow-chart of Nelder-Mead algorithm is given in Appendix B. 

Layer permeabilities in the horizontal plane and in the z-direction were set as model parameters, 
as well as their respective skin factors. Therefore, the vector of model parameters x was defined as: 


X — \^hxyl 5 kxy2 7 ... 5 kxyn , kzl , kz2 , . . . , ^zni 'S'l, *8*2, . . . , * 8 n}'. 

Layer porosities, wellbore lateral lengths, layer thicknesses and rock and fluid compressibilities 
are assumed to be constant and known. For this reason, they were not considered as model pa¬ 
rameters to be determined by the method. It was assumed that pressure and layer flow-rates were 
measured at Nme times. Thus, the objective function to be minimized was defined as: 


0{x) 


Pwf - Pwf(x) 


+ E 

1=1 


qf^ - q^ix) 


(35) 


where ||p()f /||2 is the Fuclidean norm of the Ai„e-dimensional vector that contains the observed 
pressure data at all time steps, \\p'^f — Pu,/(a :)||2 is the Fuclidean norm of the Ame-dimensional 
vector that represents the mismatch between the observed pressure and the pressure data estimated 
using vector x of model parameters, ||q '°^^||2 is the Fuclidean norm of the Ame-dimensional vector 
that contains the observed flow-rate data in layer j at all time steps and — qj{x )\\2 is the 

Fuclidean norm of the Ame-dimensional vector that represents the mismatch between the observed 
flow-rates in layer j and the flow-rate profile for layer j computed using vector x of model parameters. 

The objective function given in Fq. (35) is expected to present a unique minimum. Moreover, 
in a theoretical case where pressure and flow-rate measurements are free of noise and errors, the 
objective function definition will imply that 0{x) = 0. Field data, however, are always subject to 
noise. Therefore, even if the vector of model parameters x matches the true parameter values, some 
mismatches between computed and observed data will occur, and the objective function will be 
non-zero. Nonetheless, since the objective function defined in Fq. (35) accounts for both pressure 
and flow-rate data, it is expected to have a unique or acceptable local minimum. 
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<1 

(m3/d) 

(cP) 

Cr 

(kgf/cm^)"^ 

Co 

(kgf/cm^)"^ 

B 

(m^/m® std) 

2000 

3.1 

8.0x10-^ 

1.14x10-* 

1.0 


Table 1: Input Reservoir Parameters 


The convergence criteria was set based on the decrease of the objective function. At each itera¬ 
tion, the simplex vertex that presents the highest objective function values is replaced. Therefore, it 
is expected that, at each iteration, the maximum objective function value presented by the new sim¬ 
plex decreases. If the absolute difference between the objective function maximum value considering 
all simplex vertices is signihcantly small after two successive iterations, then further iterations will 
not improve the estimates for model parameters. In this work, we established that if the absolute 
difference between objective function maximum value after two successive iterations is smaller than 
1x10“^®, then converged was reached. 

NM was originally proposed for unconstrained optimization (Nelder and Mead, 1965). There¬ 
fore, if no restrictions are applied, the algorithm may eventually generate a simplex with negative 
values for layer permeability, which does not make physical sense. It is possible, however, to impose 
boundary constraints to the algorithm (Ghiasi et ah, 2007). A simple constraint consists of dehning 
a range that should contain the estimated parameters. This kind of linear boundary constraint is 
easily coupled to the NM work-chart depicted in appendix B by checking, at each iteration, if all 
members of the updated simplex are within the determined range. In this work, the linear boundary 
constraints were applied to NM algorithm as described by Luersen and Le Riche (2004). 

Recalling Eq. (34), layer permeabilities and skin factors were set as model parameters. For each 
layer, the permeability lower limit was set as 10 mD and the upper limit was set as 1000 mD. In 
all cases, both horizontal and vertical permeabilities for the initial simplex were randomly chosen 
inside this range. Initial skin factors, in their turn, consisted of randomly chosen values between 
-1.0 and 8.0. 


4 Results 

The techniques described in Section 3 to estimate individual layer properties were applied at a 
pair of synthetic cases to assess their accuracy. Production flow-rate, oil viscosity and rock and oil 
compressibilities were the same in both cases. The values are reported in Table 1. Layer properties 
in each case are displayed in Table 2. Effective wellbore lengths, layer porosities, compressibilities 
and the distances between wellbore laterals and a reservoir vertical boundary were assumed to be 
known. A production period of 480 hours (20 days) was considered. This test time was long enough 
for the late-time radial flow regime to develop in both cases. No boundary effects are observed, 
because the reservoir was assumed to be laterally inhnite. 

Figure 2 shows the reference pressure and pressure derivative data for both cases, while Figure 3 
exhibits the reference flow-rate data. In each case, the reference data consists of 61 distinct pressure 
points and 61 flow-rate measurements per layer. Data shown in Figures 2 and 3 were computed 
using the analytical model developed by Pan et al. (2010), as detailed in Appendix A. 


Case 

(pj 

hj 

(m) 

dZj 

(m) 

Lj 

(m) 

h 

7 

(mD) 

kzj 

(mD) 

'^ZXJ 

(mD) 

Sj 


0.16 

20 

10.0 

700 

350 

350 

350 

4.1 

A 

0.30 

20 

10.0 

800 

480 

480 

480 

2.1 


0.23 

25 

12.5 

600 

570 

570 

570 

6.1 

B 

0.24 

20 

10.0 

600 

600 

300 

424 

3.3 

0.30 

15 

7.5 

500 

500 

250 

354 

4.8 


Table 2: Layer Properties in Each Case 
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Fig. 2: Reference pressure and derivative data for cases A (left) and B (right) 


To represent a more realistic condition, an artificial noise was added to the computed pressure 
and flow-rate data. Flow-rate noise consisted of a Gaussian noise with mean zero and a standard 
deviation equal to 1% of the noise-free flow-rate in each layer. The pressure noise, in its turn, presents 
two components: one sinusoidal term and one Gaussian term. The sinusoidal term, which represents 
a tidal effect, has amplitude equal to 1% of the maximum true pressure data. The Gaussian term 
presents mean zero and standard deviation equal to 0.1% of the maximum true pressure data. 

In a real case, higher measurement errors may be observed. However, as stated in Section 3, 
the application of the RNPA and DTM requires that the constant derivative levels associated with 
early- and late-time radial flow are clearly identified. For this reason, the artificial noise was set 
with a relatively small amplitude, so that constant derivative levels may still be clearly identified. 

All three methods were applied considering the noisy pressure and flow-rate data shown in 
Figures 2 and 3. Since wellbore lengths and reservoir properties are different in each case, early- 
and late-time radial flow occur at distinct time ranges (Ozkan et ah, 1989; Odeh and Babu, 1990). 
The vertical dashed lines in Figure 2 indicate the moment when early-time radial flow ends, while 
the dotted vertical lines represent the start of late-time radial flow regime. 

Hence, the early-time logarithmic approximations presented in Eqs. (1) and (14) are valid for 
all times before the vertical dashed lines, whereas the late-time approximations given by Eqs. (7) 
and (15) are applicable for all times after the dotted vertical lines in Figure 2. In other words, it 
was assumed that early-time radial flow is simultaneously occurring in all layers until the moment 




Fig. 3: Reference flow-rate data for cases A (left) and B (right) 
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Case 

{^zx^eq 

(mD) 

Seq 

{^xy^eq 

(mD) 

Steq 

A 

464 

2.11 

520 

-7.56 

B 

413 

2.50 

605 

-7.28 


Table 3: Estimated reservoir equivalent properties 


signed with a vertical dashed lines (and analogously with late-time radial flow and the vertical 
dotted lines). 

NM, in its turn, aims to match pressure and flow-rate data for the entire test time according to 
the objective function defined in Eq. (35). For this reason, NM does not reqnire a clear identification 
of any specific flow regime. As mentioned in Section 3.4, the initial simplex is built by randomly 
sorted values within a determined range. Thus, the Nelder-Mead algorithm was performed five times 
for each case to assure that the algorithm converges to the same set of properties even if the initial 
simplex is changed. At each run, different initial guesses were used. The initial simplex parameters 
in each run consisted of randomly sorted values with uniform distribution. As mentioned in Section 
3.4, linear constraints for the NM algorithm were specified. The minimnm layer permeability value 
was set as 10 mD, whereas the maximum layer permeability value was set as 1000 mD. As for the 
layer mechanical skin factors, it was considered a range between -1.0 and 8.0. Differences between 
estimated layer properties for each trial were smaller than 1%. These results evidence that the initial 
guesses for NM did not affect the method convergence towards the unique global minimum. In this 
section, results for the first trial are reported. 

Conventional pressure transient analysis depicted in Section 2 was performed to estimate the 
reservoir equivalent parameters. Thus, the reference pressnre and pressure derivative data displayed 
in Figure 2 were used to determine the reservoir equivalent permeabilities (as indicated in Eqs (5) 
and (12)) and the reservoir mechanical and total equivalent skin factors (as mentioned in Eqs. (6) and 
(13)). The obtained results, which may be observed in Table 3, are relevant dne to the reasonableness 
check that must be performed to assess the consistency of estimated layer parameters. As described 
in in Section 3.3, if layer properties are successfully estimated, then the parameters displayed in 
Table 3 must be close to the equivalent reservoir properties obtained by RNPA and DTM. 

As explained in Section 2, at late-times, the stream lines convergence becomes relevant, inducing 
a pseudo-skin beside the mechanical damage aronnd the wellbore (Daviau et ah, 1988; Kuchuk 
et ah, 1991; Yildiz, 2003). Stream lines convergence may occur differently in each layer, making 
it unfeasible to tell apart the pseudo and mechanical skin effects dnring late-time radial flow in 
multilayer reservoirs. For this reason, the equivalent skin computed using late-time radial flow data 
represents a total equivalent skin, as defined in Eqs. (9) and (10). Estimates for the equivalent 
mechanical skin factor as defined in Eq. (4) were compnted using early-time radial flow data only. 

Table 4 displays, for each case, the objective function values considering the trne model parame¬ 
ters and the estimated layer properties obtained by the NM method. Estimates for the permeability 
in ^-direction {kzj) are also shown. It is interesting to notice that the estimated parameters yielded 
a lower objective fnnction valne than the true set of layer properties, which may be explained by 
the artificial noise added to the observed data. Estimates for kzj presented a significant error in 


Case 

0{xtrue') 0(^Xest') 

kzj (mD) 

True Est. Err. (%) 

A 

0.0496 0.0475 

350 341 -2.4 

480 461 -4.0 

570 512 -10.1 

B 

0.0324 0.0313 

300 267 -10.9 

250 214 -14.4 


Table 4: Objective function values at the final iteration of the NM method, objective function values 
considering the vector of true model parameters and estimated values of kzj in each layer 
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Case 

True 

RNPA 

DTM 

NM 1 

Est. 

Err. (%) 

Est. 

Err. (%) 

Est. 

Err. (%) 


350 

372 

6.4 

308 

-12.1 

346 

-1.3 

A 

480 

460 

-4.2 

516 

7.4 

471 

-1.8 


570 

678 

19.0 

604 

6.0 

544 

-4.6 

B 

424 

427 

0.6 

418 

-1.6 

401 

-5.4 

354 

385 

9.0 

399 

13.0 

329 

-7.0 


Table 5: Estimated vertical permeabilities (kzxj) in each layer (permeabilities in mD) 


case B. However, it is more meaningful to analyze the estimated permeabilities in the vertical plane 
{kzxj) rather than the permeabilities in ^-direction. 

Tables 5 and 6 show the estimated values for layer permeability in the za;-plane and mechanical 
skins, respectively. One should recall that the layer permeability in the za;-plane {kzxj) is not 
directly estimated by the NM method, since only layer permeabilities in the ^-direction and in the 
ccy-plane were set as model parameters. Table 5 shows the geometric mean of the values for kzj and 
kxyj computed by the NM algorithm, to make the comparison with the RNPA and DTM easier. 
Analogously to the conventional analysis, estimates for kzxj and Sj using RNPA and DTM in Table 
6 consider early-time data only, as outlined in Sections 2, 3.1 and 3.2. Skin factors computed from 
the NM algorithm, on the other hand, use data from the entire test, since this method does not rely 
on the identification of any flow regime. 

As evidenced by Table 5, the three methods yielded good estimates for layer permeabilities. 
Given that result, one should notice that the RNPA and DTM are much easier to implement than 
the NM method. Provided that early- and late-time radial flow regimes have been properly identified, 
the 2 techniques based on logarithmic approximations demand only a general worksheet software to 
be performed. The NM method, in its turn, requires not only a computational implementation of the 
analytical solution proposed by Pan et al. (2010), but also regular flow-rate measurements during 
the test. The simpler implementation of the RNPA and DTM, however, explains the higher errors 
observed in some cases (for instance, the estimated permeability of layer 3 using RNPA), since less 
data are used by these methods. Thus, the results were acceptable considering the trade-off between 
simplicity and accuracy, and that the techniques based on the logarithmic approximations use only 
part of the pressure and flow-rate history, whereas the NM algorithm uses data from the entire test. 

Even though layer permeabilities were well estimated. Table 6 shows that only the NM method 
was able to provide accurate estimates for layer mechanical skin factors. This is possibly explained 
by the fact that the two methods based on logarithmic approximations are more sensitive to data 
measurement errors (represented by the artificial noise) than the NM algorithm. As evidenced by 
Eq. (18), layer skin factors are estimated by the RNPA using the pressure and flow-rate data at 
only one point (in this work, layer skin factors were estimated considering pressure and flow-rate 
data at t = 0.006 h). DTM, in its turn, depends on data at two distinct time points, as explained 
in Section 3.2. Therefore, the two methods based on logarithmic approximations are more sensitive 
to measurement errors, despite being much simpler than the NM algorithm. 

Comparing the results in Tables 5 and 6, it is also interesting to notice errors magnitude con¬ 
sidering RNPA and DTM were similar. As detailed by Galvao and Guimaraes (2017), the need to 
repeatedly measure layer flow-rates consists of a significant operational setback of RNPA. On the 
other hand, DTM requires PLT data from only 2 distinct time steps during early-time radial flow 


Case 

True 

RNPA 

DTM 

NM 1 

Est. 

Err. (%) 

Est. 

Err. (%) 

Est. 

Err. (%) 


4.1 

2.9 

-29.2 

1.7 

-57.1 

4.1 

1.2 

A 

2.1 

1.0 

-54 7 

1.5 

-31.0 

2.1 

0.9 


6.1 

5.1 

-16.6 

4.4 

-28.2 

5.9 

-4.1 

B 

3.3 

2.0 

-39.6 

1.7 

-46.4 

3.0 

-7.4 

4.8 

3.5 

-26.2 

3.8 

-20.4 

4.3 

-9.2 


Table 6: Estimated mechanical skin factors {Sj) in each layer 
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Case 

True 

RNPA 

DTM 

NM 1 

Est. 

Err. (%) 

Est. 

Err. (%) 

Est. 

Err. (%) 


350 

399 

13.9 

396 

13.1 

350 

-0.1 

A 

480 

573 

19.4 

478 

-0.4 

482 

0.4 


570 

610 

7.0 

546 

-4.2 

578 

1.4 

B 

600 

657 

9.5 

657 

9.4 

603 

0.5 

500 

556 

11.2 

479 

-4.1 

505 

1.1 


Table 7: Estimated horizontal permeabilities {kxyj) in each layer (permeabilities in mD) 


and 2 time steps during late-time radial flow. Therefore, results from tables 5 and 6 indicate that the 
DTM is preferable to RNPA, since both yielded estimates for layer properties with similar accuracy. 

Table 7 exhibits the estimated horizontal permeabilities for all three methods. Permeabilities 
from RNPA and DTM were determined using late-time data, while NM estimates were computed 
considering pressure and flow-rates profiles during the entire test. Results show a decent accuracy, 
expect for the RNPA in case A. This is possibly related to the sinusoidal component of the artificial 
noise added to the pressure data. This oscillatory component is more relevant at late-times and, 
therefore, directly impacts the horizontal permeability estimates. Nonetheless, both the DTM and 
the NM algorithm were able to provide good estimates for layer permeabilities. Moreover, comparing 
Tables 5 and 7, the estimated layer permeabilities indicate that anisotropy effects are substantially 
more pronounced in case B than in case A, which is expected, since case A consists of a isotropic 
reservoir. 

One should notice that NM only provides one estimate for the mechanical skin in each layer, as 
it attempts to hnd a set of parameters that fit pressure and flow-rate data during the entire test. 
In contrast, the RNPA and DTM allow determination of layer total skin factors (Stj) from late¬ 
time radial flow data. Then, an alternative estimate for the mechanical skins in each layer may be 
computed from Eqs. (10) and (11). Estimates for layer total and mechanical skin factors computed 
from late-time data are reported in Table 8. 

In both cases, total skin factor for each layer was determined with decent accuracy. However, 
errors were hugely magnihed when the estimated total skin factors were applied in Eq. (10) to 
evaluate the mechanical skin in each layer. These results show that the mechanical skin computation 
based on late-time radial data is highly sensitive to errors in the determination of total skin factors. 
As dehned in Eqs. (10) and (11), estimates for both horizontal and vertical permeabilities and total 
skin factors are required to obtain the mechanical skin from late-time radial flow data. Therefore, 
computation of the mechanical skin is influenced by the accuracy of three other estimated properties, 
specially the total skin factor. Moreover, the ratio between layer wellbore length and thickness also 
magnifies the errors from total skin estimates. 

In Section 3.3, it was stated that a reasonableness check must be performed to verify if estimated 
layer parameters are consistent with the reservoir equivalent properties (Galvao and Guimaraes, 
2017; Guimaraes and Galvao, 2017), which are determined through conventional pressure analysis 
techniques. Thus, for each method, estimated layer properties were applied in Eqs. (2), (4), (8) 
and (9) and compared to the equivalent properties reported in Table 3, which were estimated 
from the pressure profiles and using Eqs. (5), (6), (12) and (13). As mentioned in Section 3.3, 



Total skin (Stj) 

Mechanical skin (Sj) I 

Case 

True 

RNPA 


DTM 

True 


RNPA 


DTM 


Est. 

Err. (%) 

Est. 

Err. (%) 

Est. 

Err. (%) 

Est. 

Err.(%) 


-7.52 

-7.51 

-0.1 

-7.56 

0.6 

4.1 

4.3 

4.7 

2.5 

-38.6 

A 

-7.73 

-7.67 

-0.7 

-7.85 

1.7 

2.1 

4.3 

104 

-3.0 

-240 


-7.16 

-7.35 

2.7 

-7.58 

5.8 

6.1 

1.6 

-74.4 

-3.8 

-162 

B 

-7.26 

-7.32 

0.8 

-7.36 

1.4 

3.3 

2.0 

-38.7 

1.1 

-67.1 

-7.05 

-7.16 

1.5 

-7.38 

4.6 

4.8 

2.3 

-50.7 

-2.9 

-160.2 


Table 8: Estimates for total skin factors {Stj) and mechanical skin factors {Sj) in each layer using 
late-time radial flow data 
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Case 

Property 

Conv. An. 

RNPA 

Estim. DifT. (%) 

Estim. 

DTM 

DifT. (%) 

Estim. 

NM 

DifT. (%) 


ij^zx ) eq 

452 

493 

8.5 

472 

3.8 

450 

-0.9 

A 

Seq 

2.00 

3.08 

54.1 

2.60 

30.3 

3.92 

96.5 


{^xy'jeq 

520 

534 

2.6 

479 

-7.9 

478 

-8.1 


Steq 

-7.56 

-7.50 

-0.8 

-7.66 

1.3 


N/A 


ij^zx ) eq 

403 

408 

1.2 

409 

1.6 

368 

-8.6 

B 

Seq 

2.37 

2.63 

11.0 

2.65 

12.0 

3.54 

49.6 

{^xy'jeq 

605 

614 

1.5 

581 

-4.0 

561 

-7.2 


Steq 

-7.28 

-7.25 

-0.4 

-7.37 

1.1 


N/A 


Table 9: Reasonableness check for both cases (permeabilities in mD) 


the reasonableness check is particularly important to assess the results obtained through RNPA 
and DTM, since these techniqnes rely on the validity of the logarithmic approximations given 
by Eqs. (14) and (15). Fnrthermore, it was also assumed that early- and late-time radial flow 
simnltaneously occur in all layers considering the time ranges defined by the vertical lines of Figure 
2. NM application, however, requires only that the analytical formulation developed by Pan et al. 
(2010) is valid. Therefore, the acceptance criteria of maximnm 10% difference proposed in Section 
3.3 was applied only at RNPA and DTM. 

Equivalent mechanical skins for RNPA and DTM were computed using the early-time estimates 
for layer skin factors because the late-time estimates presented higher errors (as displayed in Tables 
6 and 8). Fnrthermore, isolating an eqnivalent mechanical skin from the total eqnivalent skin factor 
evalnated during late-time radial flow is a complicated task. Thus, the reasonableness check for 
skin factors computed from late-time data should be performed with respect to total skin factors, 
rather than mechanical skins. Since the NM algorithm does not provide direct estimates for total 
skin factors, no reasonableness check was performed regarding the total skin factors for NM. 

The reasonableness checks are exhibited in Table 9. For both RNPA and DTM the estimated 
layer properties provided equivalent reservoir permeabilities and total skin factors are very close 
to the valnes obtained by conventional analysis, and shonld be considered valid according to the 
maximum acceptable difference established in Section 3.3. 

Nonetheless, the equivalent mechanical skin factors evaluated by RNPA and DTM presented 
significant differences compared to the valnes obtained by conventional analysis. These results are 
an immediate consequence of the poor layer mechanical skin estimates, reported in Table 6. In 
snmmary, the reasonableness checks endorse the accnracy of estimated layer permeabilities for all 
methods, but evidenced that the techniqnes based on logarithmic approximations failed to provide 
good estimates for layer mechanical skins. 

Equivalent permeabilities compnted using the individnal layer permeabilities obtained from NM 
were also close to values determined through conventional analysis. However, interestingly enough, 
the NM equivalent mechanical skins also presented relevant differences despite the fact that this 
technique provided good estimates for individual layer skin factors (as shown in Table 6). This result 
is possibly explained by the artificial noise added to the pressure data. The equivalent mechanical 
skin is determined through Eq. (13) and, thus, is also sensitive to pressnre measurement errors. 

A final accuracy assessment was performed by comparing the reference pressure and flow-rate 
data to the pressure and flow-rate profiles computed by applying the layer properties estimated 
by each method into the analytical model developed by Pan et al. (2010). Figure 4 shows the 
comparison between pressure and pressure derivative profiles for each method and the reference 
data, while Figure 5 displays the comparison between flow-rate profiles. The NM method is the 
only techniqne that aims to fit layer flow-rate data, according to the objective fnnction defined 
in Eq. (35). Therefore, Figure 5 compares the reference flow-rate data only to the flow-rate data 
obtained using the parameters determined by the NM method. In Fignre 4, RNPA and DTM 
pressure data were computed considering that layer skin factors are equal to the values portrayed 
in Table 6. That is, estimates for skin factors obtained from early-time data were used to determine 
the pressure profile. 
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Fig. 4: Comparison between estimated and reference pressure data for cases A (left) and B (right) 


A good agreement is observed between reference pressure derivative data and the pressure deriva¬ 
tive profiles computed using the estimated layer properties for all methods. Pressure derivative 
reflects the reservoir equivalent permeability (Larsen, 1982; Ehlig-Economides and Joseph, 1987; 
Galvao and Guimaraes, 2017). Therefore, the good agreement between pressure derivative curves 
observed in Figure 4 endorses that layer permeabilities were accurately estimated by all methods. 

Pressure change, in its turn, is also affected by layer skin factors. Pressure change curves com¬ 
puted using the parameters obtained by the RNPA and DTM are slightly lower than the reference 
data. The NM pressure data, on the other hand, presented a great agreement with the reference 
data. These differences, combined with the pressure derivative profiles, suggest that only the NM 
algorithm provided good estimates for layer skin factors, which is consistent with the comparison 
between estimated layer mechanical skins shown Table 6. 

Finally, Figure 5 shows an excellent match between the reference data and layer flow-rate profiles 
obtained using the parameters estimated by the NM method. This great match endorses the accuracy 
of the NM algorithm. 

As previously mentioned, the NM algorithm does not rely on the identification of any specific 
flow regime. To verify its performance in cases where early- or late-time radial flow are not clearly 
identified, the NM method was applied two other times in case A, discarding part of the reference 
pressure and flow-rate data. First, data for t > 100 h were discarded, to represent a case where only 




Fig. 5: Comparison between estimated and reference flow-rate data for cases A (left) and B (right): 
Layer 1 is represented in black, layer 2 is represented in blue and layer 3 is represented in red 
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Oix) 

k 


k^j (mD) 

k 

zxj (uiD) 

Sj 

True Est. 

True 

Est. Err. 

True 

Est. 

Err. 

True 

Est. Err. 

True 

Est. 

Err. 


350 

349 -0.4 

350 

325 

-7.3 

350 

336 -3.9 

4.1 

3.9 

-3.9 

0.0547 0.0520 

480 

480 -0.1 

480 

446 

-7.0 

480 

463 -3.6 

2.1 

2.0 

-5.6 


580 

581 1.9 

570 

483 

-15.2 

570 

530 -7.1 

6.1 

5.6 

-8.2 


Table 10: Results for case A obtained by the NM method considering only data for t < 100 h 


early-time radial flow is observed. For the second trial, pressure and flow-rate data for t < 1x10“^ 
h were discarded. Thus, only late-time radial flow may be identified. 

The results considering only data for t < 100 h may be seen in Table 10. Again, the method 
was performed five times, yielding a difference lower than 1% between the results from each run, 
and Table 10 displays the results for the first run. One should recall that some data points were 
discarded and, therefore, the objective function value computed using the true set of parameters is 
slightly different than observed in Table 4. 

Layer properties were estimated with decent accuracy. This shows that the NM method was able 
to provide good estimates for layer horizontal permeabilities despite the fact that late-time radial 
flow was not identified before 100 hours. Figure 6 compares the reference data to the pressure and 
flow-rate profiles computed using layer properties exhibited in Table 10. A great match is observed. 
Results from Table 10 and Fig. 6 evidence that layer properties could be accurately determined 
using the NM algorithm, despite the fact that late-time radial flow data were discarded. 

On the other hand, results were substantially worse when early-time data was discarded. Table 
11 exhibits the estimated layer properties at each individual method run. As evidenced in Table 
11, the method only converged in three out of the five runs, although the objective function at 
the last iteration presented similar values. Moreover, permeabilities in the z-direction and layer 
mechanical skins showed significant differences at each run. These results suggest that early-time 
data are critical to accurately determine k^j and Sj in each layer. 

Despite that, it is interesting to observe that layer horizontal permeabilities were determined 
with good precision in all runs. This is a relevant outcome, as it shows that the NM method may 
provide good estimates for kxyj even if early-time data is not available. 


5 Summary and Conclusions 


This work provides methods for determining layer permeabilities and skin in multilayer reser¬ 
voirs with multilateral wells. We extended the RNPA (rate normalized pressure analysis) (Ehlig- 




Fig. 6: Comparison between estimated and reference pressure (left) and flow-rate data (right) con¬ 
sidering only data for t < 100 h 
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Parameter 

True 

1st run 

2nd run 

3rd run 

4th run 

5th run 

Convergence 

N/A 

Conv. 

Not 

Conv. 

Conv. 

Not 

0{x) 

0.0513 

0.04937 

0.04969 

0.04963 

0.04940 

0.04970 


350 

351 

349 

349 

351 

350 

h 

f^xyj 

480 

485 

485 

485 

485 

484 

570 

582 

576 

575 

584 

578 


350 

391 

260 

260 

434 

325 

h 

i^ZXJ 

480 

569 

509 

511 

602 

449 

570 

613 

639 

565 

659 

515 


4.1 

3.5 

7.1 

7.1 

2.8 

4.7 

Sj 

2.1 

1.4 

1.4 

1.3 

1.1 

2.8 


6.1 

6.6 

4.5 

6.0 

6.1 

7.8 


Table 11: Estimated layer properties obtained by the NM method at each run considering only data 
for t > lxl0“2 h 


Economides and Joseph, 1987; Raghavan, 1989) and the DTM (delta transient method) (Galvao 
and Gnimaraes, 2017; Gnimaraes and Galvao, 2017) for mnltilayer reservoirs with mnltilateral hor¬ 
izontal wells. These methods demand only a simple spreadsheet software to be performed. The 
proposed applications of RNPA and DTM in reservoirs with multilateral horizontal wells rely on a 
clear identification of early- and late-time radial flow regimes. 

In a practical case, however, early radial flow may not occur due to wellbore storage effects and 
small reservoir thickness (Goode and Thambynayagam, 1987; Odeh and Babu, 1990). Besides, the 
test may also end before the start of late-time radial flow, which is also delayed if the wellbore 
lateral lengths are too long (Ozkan et ah, 1989; Kuchuk et ah, 1991). Thereby, the Nelder-Mead 
optimization algorithm (NM) (Nelder and Mead, 1965) was also applied to estimate layer perme¬ 
abilities and skin factors. The NM method does not rely on the identihcation of any particular flow 
regime. It requires, though, a computational implementation of the analytical model for pressure 
change developed by Pan et al. (2010) and several layer flow-rate measurements throughout the 
test. The objective function to be minimized accounts for both wellbore pressure and layer flow-rate 
profiles. 

The proposed techniques were applied on a pair of synthetic cases. To represent a more realistic 
condition, artificial noise was added to the pressure and flow-rate data. Results showed that all 
methods yielded good estimates for layer permeability and successfully identified anisotropy effects. 
However, only the NM algorithm yielded decent estimates for layer mechanical skins, while both the 
RNPA and DTM failed to provide good estimates for the individual layer mechanical skin factors. 
Late-time radial flow data provided decent estimates for layer total skins using the RNPA and DTM. 
Nonetheless, the errors are highly magnified when these total skin estimates are used to compute 
the mechanical skin factors. 

The parameters computed from RNPA and DTM presented errors with similar magnitude. It 
is important to highlight, though, that DTM requires the production logging data at only four 
distinct time steps (two at early-time and two at late-time), while the RNPA demands several flow- 
rate measurements. In this work, 6 distinct flow-rate measurements per radial flow regime (resulting 
in 12 distinct flow-rate points) were used to perform the RNPA. Thus, considering the operational 
advantage of DTM, this method is preferable to RNPA for a real field application. 

Using the NM algorithm, it was possible to obtain good estimates and also avoid generating sets 
of layer properties that presented signihcant errors despite yielding a pressure response similar to 
the reference data. Nonetheless, application of the NM method requires the layer flow-rate profiles to 
be known. Since this work used the synthetic data obtained from an analytical model, determining 
layer flow-rates at each time step was not an issue. This, however, may be an operational setback 
for the successful application of the NM algorithm in a real field case. 

Applications of the NM method were also made considering that either early- or late-time radial 
flow regime are not identified during the test. Results show that this technique is able to provide 
good estimates for layer horizontal permeabilities even under those circumstances. However, layer 
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permeabilities in the ^-direction and mechanical skin factors could not be accurately determined by 
the NM method when early-time data were discarded. 

In summary, both RNPA and DTM provided a simple way to accurately evaluate layer perme¬ 
abilities but were not able to accurately estimate layer skin factors. Since DTM presents a relevant 
operational advantage, we recommend its use for real field applications, instead of RNPA. However, 
if a detailed layer flow-rate profile is available and the radial flow regimes are not clearly observed, 
then the NM algorithm can be used. This method yielded the best estimates for layer properties 
between the three tested techniques, despite being more computationally expensive. 


Nomenclature 

B = Formation volume factor (m^/m^ std or B/STB) 
c = Compressibility ((Pa)“^ or (kgf/cm^)'^ or (psi)“^) 

dz = Smallest distance between the a wellbore lateral and a vertical boundary (m or ft) 
e = Exponent required to compute the geometric pseudo-skin in Eq. (11) 
h = Thickness (m or ft) 
k = Permeability (m^ or mD) 

L = Length of the wellbore lateral (m or ft) 

m = Constant derivative level attained during early or late radial flow 
n = Number of reservoir layers 

Nme = Number of distinct time steps when pressure and flow-rate have been measured 
O = Objective function to be minimized by the Nelder-Mead algorithm 
p = Pressure (Pa or kgf/cm^ or psi) 

Pat = Rate-normalized pressure change 
q = Flow-rate (m^/s or STB/d or m^/d) 
r = Radius (m or in or ft) 

5 = Skin factor 

Sg = Pseudo-skin that occurs due to stream lines convergence, defined in Eq. (11) 
t = Time (s or h or d) 

X = Vector of model parameters in the Nelder-Mead algorithm 
a = Coefficients required to perform the Nelder-Mead algorithm 
P = Beta transient parameter, defined in Eq. (22) 

6 = Delta transient parameter, defined in Eq. (22) 

7 = Euler constant, equal to 0.57722 
p = Viscosity (cP or Pa-s) 

(ft = Porosity 

a = Instant sources integral, defined in Eq. (A-2) 

T = Silent time integration variable 


Subscripts 

c = Contraction 
e = Expansion 
eq = Equivalent 

erf = Related to the early-time radial flow regime 
i = Simplex element index 
j = Layer index 
k = Time step index 

Irf = Related to the late-time radial flow regime 
r = Reflection 
s = Shrinkage 
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t = Total 
wf = Well face 
X = Related to x direction 
xy = Related to a:y-plane 
y = Related to y direction 
0 = Related to ^ direction 
zx = Related to za;-plane 
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A Mathematical Model 

The considered mathematical model assumes the single-phase isothermal flow of a slightly compress¬ 
ible fluid with constant viscosity and compressibility. Gravitational and wellbore storage effects are 
neglected. 

It is assumed a multilayer stratified reservoir with a multilateral horizontal well, such that each 
layer is perforated by exactly one well ramification. Since no formation crossflow is assumed, a 
multilateral well is required so that the influence of all layers is felt in pressure response. Moreover, 
it is considered that each wellbore branch is parallel to the horizontal plane. It was also assumed 
that the effective wellbore length of each branch is known. In practice, the total drilled length does 
not necessarily match the well effective length, which might be challenging to determine (Boughara 
and Reynolds, 2009; Pan et ah, 2010). 

By model hypothesis, pressure is the same in all layers, apart from the hydrostatic column. It is 
considered a transversely isotropic system, that is, permeability along the horizontal plane (kxyj) 
is constant but the vertical permeability (kzj) might be different. All computations assume that a 
consistent set of units is used. 


A.I Pressure solution for single-layer reservoirs 

The real space solution for single-layer reservoirs with horizontal wells was developed by Daviau 
et al. (1988) and Ozkan et al. (1989). They considered that the wellbore may be depicted as a 
line sink. Then, using Newman’s product and Green’s functions, pressure change is decomposed as 
the product of three terms, each one related to the transient pressure pulse propagation in a given 
direction (Rosa and Carvalho, 1989): 


Ap{x,y,z,t) = j^(j{x,y,z,t) + z^^S 
L/(pCt kzx-L/ 


(A-l) 


where kzx = y/kffykz is the geometric permeability in the zx-plane and a{x,y,z,t) is defined as: 
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dr 


In Eqs. (A-1) and (A-2), S stands for the mechanical skin factor as proposed by Hawkins (1956), 
the hydraulic diffusivity in the horizontal plane is defined as rjxy = kxy/4'y'Ct, where (j) is the reservoir 
porosity, /i denotes the fluid viscosity and ct represents the total compressibility. The diffusivity in 
^-direction {rjz) is analogously dehned. Wellbore length is denoted by L, while h stands for the 
reservoir thickness and dz represents the smallest distance between the wellbore and a vertical 
boundary. The point with coordinates {x\,y\, z\) represents the wellbore heel. 

The integral term dehned in Eq. (A-2) represents the instantaneous pressure change due to the 
oil production of all sink points along the well path, considering a constant how-rate q. Therefore, 
Eq. (A-1) shows that pressure change at the wellbore is given by the time integral of all instant 
pressure changes caused by the sink well. 

Eq. (A-1) was achieved under the assumption of uniform hux along the wellbore (Daviau et ah, 
1988; Ozkan et ah, 1989). Nonetheless, this hypothesis does not match the physical reality of how 
in horizontal wells, since it implies that pressure at the well midpoint is lower than at the well 
tips. Considering that pressure (instead of huid inhow) is uniform along the well is a more realistic 
approach (Rosa and Carvalho, 1989; Jelmert and Thompson, 1991). Pressure behavior in an inhnite 
conductivity well may be determined from Eq. (A-1) using either the equivalent pressure point 
(Ozkan et ah, 1989) or the average pressure method (Goode and Thambynayagam, 1987; Kuchuk 
et ah, 1991). 


A.2 Pressure solution for multilayer reservoirs 

The real space solution for multilayer reservoirs was achieved by Pan et al. (2010). They assumed 
that layer how-rate prohle is depicted by a step-wise constant function. Then, Eq. (A-1) may be 
combined with the superposition principle to obtain an expression for the pressure change in each 
layer at the k-th timestep: 


^Pj{x,y,z,tk) 


E 


\ kjj 4>j ctj j 


'Xj{x,y,z,t), 


for j = 1, ...,n 


(A-3) 


where n is the total number of layers and crj {x, y, z,t) denotes the integral term dehned in Eq. (A-2), 
evaluated using the properties of layer j. 

Eq. (A-3) is derived from the single-lateral solution developed by Daviau et al. (1988) and 
Ozkan et al. (1989). Thus, it assumes uniform hux at each wellbore branch. Considering an inhnite 
conductivity well, pressure change at each time step must be the same in all layers. Eurthermore, 
mass balance at the well assure that the sum of all layer how-rates must be equal to the total 
production how-rate. Thus, the following linear system may be written for a n-layer reservoir at 
each time step (Pan et ah, 2010): 


I (A-4) 

[Apj{tk) = Apj+i{tk), for j = l,...,n- 1 

At each time step, after the linear system dehned in Eq. (A-4) is solved, pressure change may 
be computed by applying Eq. (A-3) at any layer. To obtain a more accurate response, each well 
lateral may be discretized into several well segments, which increases the order of the linear system 
dehned in Eq. (A-4). 
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The solution developed by Pan et al. (2010) is, in fact, more general and also applies for multi¬ 
lateral wells in single-layer reservoirs, or multilayer systems with more than one well branch in each 
layer. However, these applications are not detailed here, since the goal of this work is to compute in¬ 
dividual layer properties assuming the wellbore-reservoir schematics portrayed in hgure 1. Although 
the analytical models for single (Daviau et ah, 1988; Ozkan et ah, 1989) and multilayer reservoirs 
(Pan et ah, 2010) account for horizontal anisotropy, permeability in the x- and y-directions were 
assumed to be equal in this work. 


B Nelder-Mead Algorithm 

The Nelder-Mead algorithm (NM) is a gradient free optimization method that aims to minimize 
a n-variables function 0{x), where x = {a;i,a; 2 , ■■■Xn} is a vector containing the n variables of the 
problem (Nelder and Mead, 1965). 

NM starts by setting a simplex composed of n -|- 1 vertices. Then, the function is evaluated 
at each simplex vertex. Then, the vertices associated with the maximum and minimum function 
values are dehned as xh and xl^ respectively. The goal at each iteration is to replace xh until a 
convergence criteria is reached (Rahmati et ah, 2013; Feng et ah, 2020). 

Depending on the objective function values, four distinct operations may be performed to update 
the simplex: reflection, expansion, contraction and shrinkage. An operation is considered successful 
if it yields a new minimum. The reflection operation is defined as (Nelder and Mead, 1965; Luersen 
and Le Riche, 2004): 


Xr = OLr)x — UtXl, (B-1) 

where x represents the centroid of all simplex points excluding xh and cx,r is a positive constant 
denoted as reflection coefficient. 

Reflection is considered successful if 0{xr) < 0{xl), that is, if the reflected vertex represents a 
new minimum. If reflection succeeds, expansion is performed. This operation is dehned as (Nelder 
and Mead, 1965; Ghiasi et ah, 2007): 


Xe = OleXr + (1 “ ae)x, (B-2) 

where ae is the expansion coefficient, which must be higher than 1. If expansion succeeds, then xh 
is replaced by x^- Otherwise, if rehection succeeds but expansion fails, xh is replaced by Xr- 

If rehection fails, then it must be determined whether the rehected vertex is a new maximum, 
that is, if 0{xr) > 0{xi), \/xi 7 ^ xh- If the rehected vertex is not a new maximum, then xh is 
replaced by Xr- Otherwise, contraction is performed. The contracted vertex is computed as (Nelder 
and Mead, 1965; Feng et ah, 2020): 


where: 


Xc — CXcXaux T (1 OCc)X^ 


(B-3) 


jxH if 0{xh) > 0{Xr) 

\ Xr otherwise 


(B-4) 


and ac is the contraction coefficient, which lies inside the interval [ 0 , 1 ]. 

If contraction succeeds, then xh is replaced by Xc- Otherwise, shrinkage is performed. This 
operation consists of multiplying all simplex vertices, expect for xl hy a coefficient between zero 
and one, the shrinkage coefficient as (Nelder and Mead, 1965; Barati, 2011): 


Xi = asXi Vx^^XL■ (B-5) 

After the simplex has been updated, the suitable convergence criteria is tested. This process is 
repeated until convergence is reached. In this work, the parameters ar, cte, etc and as, required to 
perform the simplex operations, were set as I.O, 2.0, 0.5 and 0.5, respectively. Nelder and Mead 
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(1965) and Barati (2011) applied these values and showed that they provide good convergence. The 
Nelder-Mead optimization algorithm was originally proposed for unconstrained domains (Nelder and 
Mead, 1965). However, linear constraints may be easily coupled to the algorithm (Ghiasi et ah, 2007). 
Considering that, for each variable, a range that must contain the minimum has been specihed, the 
boundary constraints are coupled to the algorithm by checking, after each operation, if the new 
simplex vertex is within the specified range. If the updated simplex vertex falls out of the bounded 
range, then it is projected to the bounded region dehned by the maximum and minimum acceptable 
values (Luersen and Le Riche, 2004). Thus, if a given simplex coordinate Xi is out of the specified 
range, the projection is made as follows: 


Xi 


if Xi > xT"" 
if Xi < xT'^ 


(B-6) 


In summary, the algorithm workflow goes as follows: 

Step 1 Dehne the initial simplex (for instance, by randomly choosing the vertices). 

Step 2 Evaluate the objective function at each simplex vertex and identify xh and x^. 

Step 3 Perform the reflection operation as defined in Eq. (B-1). If reflection succeeds, proceed to step 
4. Otherwise, proceed to step 5. 

Step 4 Perform the expansion operation, as dehned in Eq. (B-2). If expansion succeeds, update the 
simplex by replacing xu by Xe- Otherwise, update the simplex by replacing xu by Xr- Next, 
proceed to step 8. 

Step 5 If rehection yielded a new maximum, proceed to step 6. Otherwise, update the simplex by 
replacing xh by x^ and proceed to step 8. 

Step 6 Perform the contraction operation, as dehned in Eq. (B-3). If contraction succeeds, update the 
simplex by replacing xh by Xc and proceed to step 8. Otherwise, proceed to step 7. 

Step 7 Update the simplex by performing the shrinkage operation, as dehned in Eq. (B-5). Next, proceed 
to step 8. 

Step 8 Apply the convergence criteria to the updated simplex. If convergence is not reached, iterate 
again by going back to step 2. 



